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I. INTRODUCTION 

One of the important properties of the behavior of the ner- 
vous system, discovered early on that has attracted sig- 
nificant amount of attention, is the synchronization of neu- 
ronal discharges. In recent years the advent of improved ex- 
perimental techniques has provided vast amounts of new syn- 
chronization data. Concomitantly, there has been a resur- 
gence in interest and controversy concerning the functional 
relevance of synchronization. It has been established that in 
vivo cortical neurons have noisy spike trains [|| (but see in 
contrast [|^]), and that groups of neurons discharge coherently 
as found in population recordings (such as EEGs, or by arrays 
of extracellular electrodes, for a review see [Q]). These two 
facts have sparked major controversies. Firstly, does noise 
(or precise timing) in neuronal spike trains contain informa- 
tion [^,||], or is information merely due to noisy process- 
ing of an average firing rate [^[Q? Secondly, is synchro- 
nization functionally (or even statistically) significant [|l^jTl]|, 
or just an epiphenomenon [|[|? In this paper we focus on 
two different aspects of synchronization that have received 
little attention so far. Can realistic neuronal networks syn- 
chronize under the biological conditions of variable intrin- 
sic neuronal properties, and the noise-induced neuronal un- 
reliability? What kind of synchronization can be obtained, 
and what are its pertinent statistical properties? It is nec- 
essary to resolve these two questions to properly formulate 
the issues to be studied in experiment, and to analyze differ- 
ent ways of probing the experimental data. Here we focus 
our attention on the extensively-studied synchronous gamma 
oscillations in hippocampus ||l2|-|l7|]. Theoretical and com- 
putational work has shown that mutual inhibition is capable 
of synchronizing neuronal networks [18,1^. Subsequent in 
vitro experiments have convincingly established the role of 
GABA-ergic hippocampal interneurons in gamma oscillations 
0,11^. Wang & Buzsaki |^] studied the effect of current 
heterogeneity and partial connectivity on the synchronization 
of the hippocampal network. They only found strong synchro- 
nization in the gamma frequency range when the current het- 
erogeneities were small [ pOpl] ] . In strong synchronization all 
neurons in a local circuit spike within a short interval of each 
other This suggests that strong synchronization can only be 
obtained when the intrinsic properties of the neurons are not 
too different. According to Ref. |^] this would mean a less 
than 10% difference in current drive, or average firing rate. It 
has been hard to pinpoint the amount of variability in intrinsic 
properties in the in vitro and in vivo preparations of differ- 
ent brain areas. It is however not unreasonable to assume the 
presence of more than 10% variability in these preparations. 
Strong synchronization is also not robust against noise [^]. 
It would therefore seem unlikely for strong synchronization 
to be present in hippocampus under physiological conditions. 
Indeed, here we show that stochastic weak synchronization 
(SWS) is more prevalent in parameter space, and is also ro- 
bust against neuronal heterogeneities and synaptic noise. We 
conjecture that as a consequence it is much more likely to oc- 
cur in neuronal systems. 



In SWS, neurons spike within a short interval from each 
other, but not necessarily at each period [23 ,2^. The syn- 
chronization is called stochastic, because the particular cycle 
in which the neurons fire is random. This makes the properties 
of this state different from the well known cluster states stud- 
ied by previous authors [^^-^sj]. There each neuron always 
fires at the same cycle with the same cluster Both strong and 
stochastic weak synchronization yield periodic population os- 
cillations. The difference can then only be ascertained using 
multi-unit recordings. 

We use cross correlation analysis to show that noise and 
heterogeneity affect the synchronization properties of our net- 
work in very different ways. Large enough noise and het- 
erogeneity will, however, stop strongly synchronized oscilla- 
tions. We demonstrate that neither adding a periodic drive nor 
increasing synaptic coupling can significantly increase robust- 
ness of strong synchronization. Finally we determine for what 
parameters robust self-induced 40Hz synchronous oscillations 
can be obtained. 



II. METHODS 

A. Single neuron model 

Our aim here is to establish physiological criteria for ro- 
bust synchronization in the gamma frequency range. The use 
of a biophysically realistic model is therefore of pivotal im- 
portance. At the same time it is also important to balance 
the amount of complexity versus practical simplicity p^]. 
We have therefore not attempted to use the latest available 
data to construct a detailed multi-compartmental model. The 
computer requirements to sample the full relevant parameter 
space, and perform our type of analysis, would be extremely 
demanding even using very fast computers. It has been shown, 
nonetheless, that one and two compartmental models can ac- 
curately generate spike trains of the right shape and frequency 
[|30|-p4]. Multi-compartmental models may be necessary to 
assess the synaptic integration of inputs located on different 
parts of the dendritic tree. This is currently an intensely stud- 
ied area in electrophysiology [|3^- |3^] . Here we study a model 
previously introduced by others [Effl. The model has been 
shown to reproduce the salient features of the dynamics of hip- 
pocampal interneurons. The neurons are modeled as a single 
compartment with Hodgkin-Huxley type sodium and potas- 
sium channels. In this work all the neurons are connected to 
all others and themselves (ALL to ALL connectivity) via in- 
hibitory GABA^-synapses. The equation for the membrane 
potential of a neuron is (the index i of the neuron is omitted) 
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Here we use: the leak current J^, = gL{V — E^), the sodium 
current iNa = gNam^h{y — Ej^a), the potassium current: 
Ik = gKn^{y ~ Ek), and the synaptic current: Igyn = 
9syns{V — Esyn)- The Gaussian noise is denoted as ^ (see 
below), and / is the tonic drive. The channel kinetics are given 
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in terms of m, n, and h. They satisfy the following first order 
kinetics: 



dx 
'dt 



0(a2;(l - x) - Pxx). 



(2) 



Here x labels the different kinetic variables m, n, and h, 
and = 5 is a dimensionless time-scale that can be used to 
tune the temperature-dependent speed with which the chan- 
nels open or close. The rate constants are poj], 



-0.1(1/ + 35) 
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exp(-0.1(y + 35)) - 1' 
4exp(-(y + 60)/18), 
0.07exp(-(F + 58)/20), 
1 

exp(-0.1(y + 28)) + l' 

-0.01(T/ + 34) 
exp(-0.1(T/ + 34)) - 1' 
0.125 exp(-(y + 44)/80). 



We make the approximation that m follows the asymp- 
totic value moo{V{t))) = a„i/{am + (3m) instantaneously. 
The synap tic gating variable s obeys the following equation 



ds 

— ^aF{Vp){l-s)-f3s, 



(3) 



with a = 12tos-i,/?= l/Tsyn,F{Vp) = l/(exp(~t/p/2) + 
1), and Vp is the presynaptic potential. The function F{Vp) is 
chosen such that when the presynaptic neuron fires, > 0, 
the synaptic channel opens. The decay time of the postsy- 
naptic hyperpolarization is chosen as Tsyn — ^/P — 10 ms 
(or 20 ms in some instances). We use a reversal potential of 
Esyn — —75 mV for the inhibitory (GABA^) synapses p^]. 
The standard set of values for the conductances used in this 



work is gNa = 35, gK = 9, gL — 0.1, and gs 



0.1 



(in mS/cm ), and we have taken En a = 55 mV, Ek = 
—90 mV, and El = — 65 mV . The membrane capacitance 
is Cm = ^nF/crn^. Unless stated otherwise we will use 
the standard set of parameters listed above. When no current 
value is specified we use / = 1 jiA/cm?. The network will 
then spike at approximately 39 Hz. 

We chose the initial values for the membrane potential at the 
start of the simulations uniformly random between — 70 and 
—50 mV . The kinetic variables to, n, h, and s are set to 
their asymptotic stationary values corresponding to that start- 
ing value of the membrane potential. 

The resulting equations with noise are integrated using an 
adapted second order Runge-Kutta method [p^, with time 
step dt =0.01 ms. The accuracy of this integration method 
was checked for the dynamical equations without noise (D = 
0) by varying dt and comparing the result to the one obtained 
with the standard 4th order Runge-Kutta method with a 
time-step dt of 0.05 ms. 

We normalize all quantities by the surface area of the neu- 
ron. This leads to the following system of units: the mem- 
brane potential V in mV, time t in ms, firing rate / in Hz, 



membrane capacitance Cm in nF/cm?, conductance g^ in 
mS/crn?, voltage noise ^ in mV/ms, strength of neuroelec- 
tric noise D in mV^/ms, the rate constants ax and Px in 
ms~^, and the current / in fiA/cm?. The kinetic variables 
TO, 71, h, s, and the time-scale <j) are dimensionless. Results in 
our paper are expressed in this system of units. 



B. Heterogeneity and synaptic noise 

We have included heterogeneity in the applied current. For 
each run we draw the applied current for each neuron from a 
uniform distribution. The average of the current distribution is 
/ and the variance is aj. The current heterogeneity represents 
the variation in the intrinsic properties of the neurons in the 
hippocampus. Experimental measurements of quantities like 
the input resistance i?i„, the membrane time-scale, the spon- 
taneous spiking rate, the shape of the somatic action potential 
(amplitude, width, rise and fall time), and the afterhyperpo- 
larization, show considerable variance [p4[-p3[]. It is hard to 
determine how much of the variance is due to measurement 
errors, and how much is actually due to intrinsic neuronal vari- 
ability. Here we assume the main effect of the variability is to 
change the intrinsic frequency of the neurons (which can be 
varied using the current drive in our model). Another source 
of heterogeneity in in vitro experiments is the glutamate pres- 
sure ejection method [|l^. It can lead to an inhomogeneous 
activation of metabotropic glutamate receptors, and thus to a 
variable current. In this paper we will consider ai as a free 
parameter 

At least three sources of noise can be identified [Q|: ran- 
dom inhibitory postsynaptic potentials (IPSP) and excitatory 
postsynaptic potentials (EPSP), stochasticity of the synap- 
tic transmission, and the stochasticity of the channel dynam- 
ics. Here we assume that the variability in the neuronal dis- 
charge is mainly due to synaptic noise [^. We have com- 
pared the effects of Poisson distributed spike trains of EPSPs 
and IPSPs to that of a Gaussian noise current on interspike 
interval (ISI) variability. Poisson and Gaussian noises, do 
not yield identical results. The statistics obtained from both 
models, however, are similar in the parameter regime stud- 
ied [p^. For the purpose of our studies we consider that 
Poisson and Gaussian distributions are two alternate ways of 
producing noisy spike trains with particular statistics. There- 
fore, the synaptic noise is only implemented as a Gaussian 
distributed, white noise current in neuron i, with {£,i{t)) = 0, 
and — 2D5{t — t')dij. The noise currents in 

different neurons are assumed independent. 



C. Measured quantities 

From our simulations we obtain the time trace for the mem- 
brane potential Vi (<) of each neuron. We determine the spike- 
trace Xi from Vi as follows: Xi{t) — 1 when Vi{t) crosses 
OmV (i.e. V{t^) < < V(t+)), and it is zero elsewhere. 
From Xi we obtain X{t) = Xi{t). X is proportional to 
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the instantaneous firing rate of the network. We also calculate 
the correlations function k pofl : 



(4) 



This function measures the amount of strong synchronization, 
and depends on the bin size r of the time discretization 
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J{XdnT)){X,{nT)) 
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We use T = 200 dt = 2 ms for oscillations in the gamma- 
frequency range, or T/10 for periodic drives with period T. 

We also evaluate other measures that yield further detailed 
quantitative characterization of the network behavior We cal- 
culate the time autocorrelation function: 



9x{t) 



{x{t)xm - (x) 



(xy) 



and the cross correlation function: 

{xit)ym~{xit)){ym 



9xy{t) 



(6) 



(7) 



Here x and y can be any of the variables Xi, Vi, and X, and 
is a shorthand notation for the time-average. 

We also consider the more conventional interspike inter- 
val histogram (ISIH) [p7|], averaged over all network neurons. 
From the ISIH one can obtain two statistics: the average ISI, 
Tjsi, and the standard deviation of the ISI, (Tjsi- The ratio 
aisi/risi is known as the coefficient of variation (CV). The 
average firing rate is / = 1 /tjsi, and the population standard 
deviation of / is cr/. 



(8) 



where fj — 1/ rjgj is the average firing rate of the jth neuron. 
In addition we plot rastergrams, with the action potential of 
each neuron plotted as a filled circle, with the y-coordinate 
given by the neuron index and the x-coordinate by the spiking 
time. 

To analyze the stochastic weak synchronization network 
dynamics we need to apply a different method. The popu- 
lation period t„ is different from the population averaged ISI, 
and to estimate it we proceed as follows. First we determine 
the firing rate X{t) as before with 1 ms bins. In the stochastic 
weak synchronization state X{t) will consist of a number of 
approximately equidistant peaks of finite width (see Fig. [lljf). 
We use the position of the first maximum of the Fourier 
transform at nonzero frequency as an estimate T for the pe- 
riod T„. We calculated the weight {X{t)), the average posi- 



T of 



tion tl = {tX{t)), and the width = ^{t^X(t)) - (4)2 
of the ith peak. The time average is taken over a range 
[-0.35r,0.35r] about the estimated position t^'^ 



the peak. We calculated the number of spikes that fall out- 
side this region. If the average number of missed spikes is 
more than one per cycle we reject the cluster state. The cycle 
length (time between two consecutive cluster firings) is de- 
fined as t'^ ^t\ — t]r^ ■ We determine the average cluster size 
= {Ni), its CV{Nc) = ^/{{NlY) -~N^/Nc, the average 
cycle length r„ = (r^), its Cl/(r„) = ^/MW^/rn, 
and the average width CTc = (Cc)- Here the average (•) is 
given by the sum over all cycles in the run (after discarding a 
transient). We characterize the strength of the synchronization 
using a modified kw and CVw- In the SWS state the ISIH 
has multiple peaks. The CV of the ISI receives contributions 
from the variance within each peak, but also of the variance 
between the multiple peaks. We are only interested in the for- 
mer, and the conventional CV is thus an overestimate. Instead 
we use CVw = <^c/ Tn which is related to the average width of 
one peak in the ISIH. The coherence k measures the number 
of coincident spikes between two spike trains. Consider two 
neurons that do not spike at each cycle, but when they both do, 
the spikes are coincident (that is in the same bin). If the prob- 
ability of spiking in a cycle is p ~ Nc/N, and both neurons 
fire statistically independent, we obtain k — p. These neurons 
can be considered synchronous and we want kw = 1- We 
therefore normalize k by p. 

There is a subtlety in the calculation of the average firing 
rate. In the deterministic noiseless case one ISI is enough to 
determine the average value (after discarding the transient). 
[Note that counting the number of spikes in a fixed interval 
is not an efficient way to determine the exact firing rate.] In 
the presence of noise, however, you need at least 10 ISIs to 
accurately determine the average. In networks with large cur- 
rent heterogeneities there are neurons with high and very low 
firing rates (Fig. [ll]). The average ISI for the low firing rate 
is less accurate than for the high firing rate neurons in the net- 
work. However, it carries equal weight in the conventional 
average tjsi = J2j '''isi- We have therefore used a weighted 
average tjsi — ''^j'''isi/ '^j (here rij is the number of 
intervals over which rjgj is calculated), and the approximate 
identity Nc/rn ~ Ns/tjsi can be used as a check. Ng is 
the number of active neurons, defined as the neurons that have 
more than two ISIs after the transient. 



III. RESULTS 

A. Non robustness of strong synchronization 

In this section we describe the results of our simulations 
for a network of = 100 interneurons, connected all to 
all, with either synaptic noise (SN), or current heterogeneities 
(CH). In Fig. ll we plot coherence parameter k (defined in 
Eq. (^) versus the strength of the synaptic noise D, and ver- 
sus the standard deviation of the current heterogeneities aj. 
We find that strong synchronization is lost for approximately 
D > 0.10 mV^/ms and ctj > 0.1 iiA/cm? (with the stan- 
dard set of parameters listed in Methods). The mechanism by 
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which strong synchronization is lost, however, is different in 
the CH case compared to the mechanism with SN. This dif- 
ference shows up only if we studies the whole state of the net- 
work using cross correlation functions, instead of the average 
quantities shown in Fig. Next we compare these two mech- 
anisms. Wang&Buzsaki (WB) [ ^o[ ] have already analyzed the 
case with current heterogeneity. We have reproduced part of 
their work, and we will refer to their corresponding figures. In 
both CH and SN cases the neuronal firing rate decreases when 
the network desynchronizes. We have plotted the time-trace 
of the synaptic drive s{t) in Figs. ^ and d. The phasic part de- 
creases, and the tonic part of s{t) increases with increasing D 
and aj. The increased tonic part is responsible for the lower 
average firing rate. The firing rate of the CH neurons satu- 
rates (when averaged over enough realizations of the current 
heterogeneities^]), whereas for the SN it increases steadily as 
a function of D for large values of D. This is because the sin- 
gle neuron firing rate increases with D [p^, but tonic inhibi- 
tion saturates to its highest value in the asynchronous network. 
The dispersion cr/ (see Methods) with CH is larger than the 
one in SN (not shown). In SN all the neurons have identical 
intrinsic properties, and the expectation value for the average 
frequency of each neuron is the same. The dispersion ct/ in 
this case represents the fluctuations in the average ISI due to 
the finite averaging time. With CH the neurons have different 
intrinsic frequencies, and the dispersion a / increases with cr/ 
(and does not go to zero after a long averaging time, see WB 
Fig. 5B). In Fig ^ and II we compare the correlation func- 
tions for the SN and CH case, respectively. In (a) we have the 
strongly synchronous network, in (c) the asynchronous net- 
work, and in (b) a transition state. The difference between SN 
and CH becomes clear when one considers the cross correla- 
tion functions. With CH the number of pairs that are phase 
locked drops gradually (see WB Fig. 8E). The pairs that are 
phase locked, are tightly phase locked (Fig. ^b 1,4,5, 6), and 
there is no dispersion in the cross correlations, only a rela- 
tive phase. Even in the asynchronous state the autocorrelation 
function for a single neuron is sharp, i.e. the neuron fires 
regularly with a fixed frequency (Fig. ^cl and c 13). The pop- 
ulation average of gx^, however, is disordered (Fig. ||lcl2), 
since each of the neurons has a different firing rate. In the SN 
case there is already dispersion due to the noise-induced jitter 
in the spike time, in the autocorrelation Fig. |^IIal 1, and in the 
cross correlations (Fig. ^Ia2-al0). The dispersion increases 
gradually with D. The difference between the CH and SN 
cases is also evident in the distribution of k values for each 
pair in the network (Fig. ||, WB Fig. 8E). For SN there is a 
well defined peak, with the average shifting to lower values 
as D increases, (Fig. pp-c), whereas for CH there is a broad 



distribution for small aj (Fig. ^), a peak at low values of k 
combined with a broad distribution for moderate values of cr/ 
(Fig. ||e). For higher values of ai the network is in an asyn- 
chronous regime, and only the peak for low k values is present 
(Fig. 1^). We have compared the ISIH for a network neuron to 
the ISIH of an isolated neuron (Figs. 0c-e), and also the values 
of Tjsi and aisi (Figs. ^ and b). The CV of the network neu- 
ron is higher than the CV of an isolated neuron which in turn 
is higher than the CV of an isolated neuron with autosynaptic 
feedback. The inhibitory coupling in the network increases 
the effect of the noise compared to uncoupled neurons: the 
jitter in the spike times reduces the phasic component of s{t) 
(Fig. ^). This effect does not take place in a neuron with au- 
tosynaptic feedback: the size of the phasic component does 
not decrease with D, only the timing deteriorates. 



B. Weak enhancement of robustness due to resonance effect 

In this section we drive the tonically active neurons with an 
external periodic drive. This drive may represent the effects 
of putative pacemaker neurons, similar to the ones that were 
recently found in the striate cortex of cats [Q. There is no 
compelling evidence for having the 40Hz pacemaker neurons 
projecting to hippocampal interneurons. Here we model the 
pacemaker as an excitatory synapse driven by a periodic pulse 
train, and investigate its effects on the robustness of synchro- 
nization. When a single neuron is driven by a periodic drive it 
will entrain, or phase lock, when the drive frequency is close 
to the natural frequency (or in some cases close to a ratio- 
nal fraction) [^]. Here we vary the natural frequency via the 
current. In Fig. 5c we show the coefficient of variation (CV) 
of an isolated neuron. When the neuron is entrained the CV 
drops to zero. The ISIH then consists of a single peak. The 
entrainment occurs for a range of current values, / = 0.87- 
1.0. Since the firing rate is constant, the f — I has a flat step 
(not shown). Outside the range of entrainment the ISIH has 
more structure (Fig. ^ and f). When there is noise present in 
the neuron, the CV will increase. For weak noise the CV in 
the entrainment regime will still be lower compared to the C V 
when the neuron is not entrained. For the network the CV is 
also lower in the entrainment regime (Fig. Ijb), though the CV 
increases faster with the noise strength compared to the CV 
of the isolated neuron. The synchronization of the network, 
measured by k (Fig. ^), is significantly enhanced in the re- 
gion of entrainment (for D = 0.004-0.04). The enhanced 
synchronization disappears for higher values D > 0.2. 



^This assumes 



f{l)dl ^ fih.) 



lav -ai^/Vi 



C. Effect of synaptic coupling strength on robustness of strong 
synchronization 

We have also studied the effect of varying the synaptic cou- 
pling strength gsyn- For large enough D the network will 
be asynchronous. We find that the network frequency in that 
case decreases with increasing values of gsyn (Fig. For an 
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asynchronous network the synaptic drive has a constant tonic 
hyperpolarizing conductance, decreasing the firing rate. The 
stronger the coupHng the larger the decrease. The synchro- 
nization measured by the parameter k displays a different be- 
havior. In Fig. we plot the k versus ggyn for one specific 
value D = 0.02 (and / chosen such that the firing rate is ap- 
proximately 39 Hz). It is interesting to note that stronger cou- 
pling does not necessarily mean a higher value of k. The co- 
herence K has a local maximum for gsyn = 0.1, for higher val- 
ues of gsyn, K decreases (see WB Fig. 12B). For gsyn > 0.3, 
K starts increasing again. We have studied the underlying dy- 
namics of this non monotonous behavior. In Fig. ^-f we plot 
the ISIH for different values of gsj^„. For gsyn > 0.2 one finds 
more than one peak. In the rastergrams (Figs. |^a and (b)), one 
can see that the dynamics corresponds to a population that has 
a well defined frequency, but individual neurons sometimes 
miss, or skip, a period. Despite this small asynchrony when 
the neuron fires, it does so in synchrony with the others. As a 
consequence the rastergrams looks much more ordered com- 
pared to the one for gsyn = 0.1 at the same noise strength 
D = 0.002. 



D. Larger gsyn leads to robust stochastic weak synchronization 

In this subsection we discuss the robust 40 Hz rhythms 
found for higher gsyn values. We have doubled the synaptic 
decay constant to Tgyn — 20. Here we will use the modified 
Kw and CVw as mentioned in the Methods section. 

In Fig. ^ we vary gsyn from 0.05 to 2.5 with a spacing of 
0.05. The neuron number is kept equal to = 100, and we 
use / = 2.0, and aj = 0. For D = 0.0 and aj ^ the 
network is in a strongly synchronized state, with the network 
frequency /„ the same as the single neuron firing rate /. The 
frequency is exactly the same as the one for a single neuron 
with autosynaptic feedback, as one would expect. This coher- 
ent state can be arrived at from many different random initial 
conditions. 

For weak noise, D = 0.008, the network stays in a strongly 
synchronized state for gsj,„ < 0.25. For higher (7sj,„ skipping 
starts to occur, the fractional cluster size decreases from val- 
ues close to one to values below one-half at gsyn — 1-2. At 
that point the network is in a real (albeit stochastic) cluster 
state, on average the neuron only fires once every two cycles. 
We will refer to all states for which certain active neurons do 
not fire at each cycle as a stochastic weak synchronized (S WS) 
network. The network frequency, /„, and the single neuron 
firing rate / both decrease with increasing gsyn- When the 
network settles in the SWS state / starts to differ consider- 
ably from its value at the D = state. The strength of syn- 
chronization increases with gsyn, that is CVw decreases and 
Kw increases. For values gsyn > 2.0, CVw and kw slowly 
saturate. 

For stronger noise D = 0.04 and D ~ 0.20 the network 
is asynchronous for low values of gsyn- We have therefore 
excluded these points based on the criterium discussed in the 
Methods section. The network frequency starts out at a higher 



value, and the neurons fire at a lower rate compared to the 
D = 0.008 case. The strength of synchronization, kw and 
CVw, is reduced compared to the one for D = 0.008, but still 
increases with gsyn- Note that all the neurons in the network 
still have a nonzero firing rate. It is thus possible to obtain 
weakly synchronized oscillations in a network consisting of 
100 neurons in the frequency range between 20 and 40 Hz. 

We find that noise is necessary to obtain SWS. We have 
studied SWS in the presence of weak current heterogeneities, 
say for aj — 0.02. Without noise (D = 0) the network is 
in an strongly synchronized state, and k displays a maximum 
as a function of gsyn (WB Figl2B). One also clearly notices 
the effect of suppression [|T]|: for larger gsyn the inhibition of 
faster spiking neurons stops the firing of neurons driven by a 
smaller current. As a result the total number Ns of active neu- 
rons gradually drops (Fig. For a small amount of noise, 
D = 0.008, the situation changes dramatically. A SWS state 
is obtained, and all neurons remain active (Ns/N = 1), while 
K saturates for gsyn > 0.5, and the value of k for gsyn > 0.8 
is even higher than without noise. Thus noise may actually 
improve the coincidence. Of course noise does increase the 
width, CVw, of the peaks in the instantaneous firing rate. The 
single neuron firing rate decreases significantly compared to 
that for L* = 0, whereas the network frequency is only weakly 
affected. 

We have performed numerical simulations for system sizes 
N = 10, 20, 50, 100, 200, and 1000. We have used the fol- 
lowing parameters values: / = 5.0, D = 0.2, gsyn = 1, 
Tsyn = 20, and ct/ = 0.1. The network frequency increases 
with system size, whereas the firing rate stays approximately 
constant with a dip around N ~ 50. The measures for co- 
herence, Kw and CVw, are also only weakly dependent on 
system size. 

The cycle-to-cycle fluctuations in cluster size vary approxi- 
mately as ViV (Fig. [Io|g). The strength of the inhibition is de- 
termined by the number of neurons that fired in the previous 
cluster, and in turn it determines at what time the first neurons 
become disinhibited. One therefore expects cluster size fluc- 
tuations and cycle length fluctuations to be intimately related. 
Indeed, the standard deviation of cycle length varies as 1/Vn 
with N the number of neurons (Fig. [lo|h). This means that 
larger networks are better at generating a precise cycle length, 
whereas size does not matter as much for the coincidence of 
spikes measured by k.w and CVw - 

In each simulation we randomly draw a set of driving cur- 
rents Ij for each neuron j from a uniform probability distri- 
bution. The results one obtains may critically depend on the 
particular realization of driving currents. One expects that for 
larger systems this is less of a problem. The population distri- 
bution of / is more likely to approach the original ensemble 
distribution of currents for a given neuron. Here we have stud- 
ied the range of values for the measured quantities (/„, /, and 
so on) for ten different realizations. We find that for most 
quantities (for these parameter values) the range of values de- 
creases with N, and for N > 500 one realization will give a 
result close to the expectation value. 

We now vary ai and D for the following fixed parameter 
set N = 1000, gsyn = 2, Tsyn = 20, and / = 3.5. For 
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D = and ct/ = the network is strongly synchronized at 
20 Hz. The instantaneous firing rate consists of a sequence 
of regularly spaced delta functions (Fig. pl]e), the ISIH has a 
single delta peak at 50 ms (Fig. [lip), and all neurons spike 
at the same frequency (Fig. pl]a). Increasing D increases the 
network frequency, but decreases the single neuron firing rate 
(Fig. The population activity is still periodic (Fig. pl]f), 
but the peaks have a finite width (as well as the ISIH), and the 
ISIH becomes multimodal. This process continues with ISIH 
spreading out more and more, with the CVw increasing, and 
Kw decreasing. 

As mentioned before we need some noise to generate an 
SWS state. Here we use D — 0.2, while at the same time 
varying ct/. For finite aj there is still a coherent population 
activity (Fig. ITllg ), despite the fact that neurons have different 
firing rates (Fig. \n\:). Increasing aj will reduce coherence, 
Kw decreases and CVw increases. At the same time both /„ 
and / increase (Fig. p^). This is different from the effect of 
increasing D. Higher aj leads to suppression, with fast spik- 
ing neurons preventing slower ones from firing, and as a result 
part of the inhibition disappears, while further increasing the 
firing rate and its average (calculated from the active neurons). 
On the other hand noise increases the tonic inhibition for each 
neuron, and thus leads to a reduced firing rate. Also the pro- 
gression of the asynchronous state is different. The first peak 
in the ISIH becomes broader, and the higher order ones have 
a reduced prominence (Fig. 111). For increasing D the peaks 
just wash out. 



E. Comparison to externally induced stochastic weak 
synchronization 

We find that robustness can be enhanced by getting a net- 
work in a SWS state. The single neuron discharge in an SWS 
state is very similar to that obtained in Stochastic Resonance 
[ |5C| ] using a subthreshold external drive. In this subsection 
we compare the previous self-induced SWS state to the one 
induced by an external drive. We drive the system by a si- 
nusoidal current of amplitude 1.2 (admittedly large). The av- 
erage value of the driving current is zero. For weak noise 
the spikes of a single neuron are spaced many cycles apart 
(Fig. pj|a). With increasing noise the ISIH starts to look more 
like the ones in Fig. [ll|b. The population activity also is pe- 
riodic (Fig. pj|e). The distribution of tjsi in Fig. pj| a looks 
similar to the one in Fig. [llja. For high values of D the neu- 
ron can spike more than once during one cycle, as a result 
the number of coincident spikes is reduced. The addition of 
current heterogeneity does not seem to affect the network be- 
havior (as long as all the neurons are still below threshold). 
The different ISIH in Fig. pj] ! look very similar. However the 
corresponding range of tjsi values in Fig. |ll]c does increase 
with a J . 



IV. DISCUSSION 

Previous authors have recognized that strong synchroniza- 
tion is only moderately robust against neuronal heterogene- 
ity |20 2l||. We have previously shown that the same holds 
including synaptic noise [p2[|. The basic premise of syn- 
chronization by mutual inhibition is almost trivial, since the 
network consists of intrinsically periodically spiking neurons. 
Their output produces a periodic synaptic drive, which in turn 
is fed back into the network. Inhibition thus allows a phase 
lock at zero phase with this drive. Heterogeneity and noise 
reduces the phasic, and increases the tonic part of the synaptic 
drive, leading to a reduction in synchronization, and eventu- 
ally leading to an asynchronous state (Fig. [l|c,d). The syn- 
chronization behavior of networks of physiological realistic 
neurons, however, is by no means fully understood. In this 
work we showed that the loss of synchronization proceeds via 
different mechanisms in the presence of synaptic noise com- 
pared to the presence of current heterogeneity. This is evident 
from the cross correlations shown in Figs, gand |^. We also 
found that the noise-induced precision loss in the uncoupled 
neuron is exacerbated by the inhibitory coupling. All of these 
could seem obvious based on previous work on heterogene- 
ity [ ^olpl] ]. However, its consequences for real life biolog- 
ical networks had not been fully appreciated. Our results, 
combined with previous results, show that there is a prob- 
lem with strong synchronization by mutual inhibition, since 
it is unlikely to occur in in vitro or in vivo systems. [There 
are exceptions such as for example the pacemaker nucleus 
in electric fish [|T]|, where the neurons are coupled via gap 
junctions.] The aim of this paper was to treat this problem: 
how can one obtain robust synchronization in the presence of 
synaptic noise and neuronal heterogeneity? Our results are 
twofold. First, methods to increase the robustness of strong 
synchronization have been ineffective. Second, we showed 
that robust stochastic weak synchronization (SWS) can be ob- 
tained for biophysically realistic parameter values. SWS is 
consistent with previous experimental data. In what follows 
we discuss these two important results in more detail. 

We believe that strong synchronization is not robust 
enough. To make sure that we do not prematurely discard 
strong synchronization by mutual inhibition we have to make 
an effort to increase robustness. In this paper we discussed 
two simple methods to increase robustness of strong synchro- 
nization. One method was to increase the synaptic coupling 
gsyn- Inhibition is responsible for synchronization. It is then 
quite natural to expect that increasing the strength of inhibi- 
tion increases robustness. The fact that this does not happen 
is surprising. For current heterogeneity this is in part due to 
suppression pl] , p0| ]. We have studied this effect for synaptic 
noise in more detail. We found that neurons skip periods for 
higher values of gsyn (see Fig. In other words the strongly 
synchronized state becomes unstable, and a weakly synchro- 
nized state emerges. This weakly synchronized state looked 
more coherent (Fig. 0a,b), and provided the impetus to further 
study the robustness of the SWS states. 

We also added a periodic drive to the neuron. Recent exper- 
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imental work shows that the CV of neurons on an entrainment 
step is reduced compared to the CV outside the step [^. A 
clear physiological correlate in hippocampus of this drive is 
lacking at present. Since we try to reject our conjecture this 
lack of physiological realism is not a problem. We found a 
moderate increase in robustness (see Fig. H). The periodic 
drive is not as effective as one would have intuited, how- 
ever In fact the inhibitory connections reduce the increase 
in robustness compared to the increase in the single uncou- 
pled neuron (see Fig. If we add a subthreshold periodic 
drive with noise to a quiescent neuron we obtain weak syn- 
chronization. This is known as Stochastic Resonance in ex- 
citable systems ||50|]. To summarize: our attempts to signif- 
icantly increase robustness of strong synchronization failed. 
Instead we found weak synchronization. Thus it is easier to 
find weak synchronization in parameter space than it is to find 
strong synchronization. 

If one accepts the fact, however, that strong synchronization 
is not robust against noise and heterogeneity, and that periodic 
population oscillations are found in experiments, then one has 
to carefully consider the possible relevance of weak synchro- 
nization. Weak synchronization as well as strong synchroniza- 
tion lead to a periodic population discharge, and specifically 
to an inhibitory synaptic drive indistinguishable from the one 
found in pyramidal neurons in [[l4[|. Moreover, the clusters 
that form in stochastic weak synchronization bear a resem- 
blance to the neuronal assemblies found in some experiments 
1 52], and that are thought to play a role in putative binding 
1 53]. The question then is, is SWS more robust, and can it 
be found for the gamma frequency range for biophysically re- 
alistic parameters? What is needed is a higher total synaptic 
conductance, and noise. The necessary amount of noise is 
very small, D > 0.004 is sufficient. The noise prevents the 
occurrence of suppression (Fig. H). In suppression the faster 
neurons prevent the slower ones from firing. This reduces the 
inhibition of the faster ones, and allows them to fire at differ- 
ent frequencies and at random relative phases. Suppression is 
thus detrimental to synchronization. 

We obtained SWS for different system sizes (we studied 
networks from 10 to 1000 neurons). The coincidence proper- 
ties {kw, and CVw) did not vary much with size. The tempo- 
ral precision of the population oscillation, however, increases 
approximately as 1/ ^/N (Fig. [l0|h). Large networks can thus 
produce precise pacemaker rhythms. In addition the statistical 
quantities in small networks show more variation with differ- 
ent realizations of the current drive. It is of considerable inter- 
est to understand why weak synchronization is so much more 
robust and prevalent compared to strong synchronization. In 
strong synchronization one requires an equal firing rate for 
each neuron, and coincident spikes, while weak synchroniza- 
tion requires only coincident spikes. By definition, then, weak 
synchronization is easier to generate. In fact strong synchro- 
nization is intrinsically less robust, for it requires phase lock- 
ing between neuron pairs. This is only possible (depending on 
intrinsic properties) for a small difference in driving currents. 
There is a price to pay, there will be a small phase difference 
between the firings of each neuron. Pairs with a large phase 
difference are less stable against the influence of noise. It is 



instructive to consider this problem as the nonlinear dynam- 
ics of a single neuron driven by a periodic drive. The neuron 
can be entrained on different n:m steps. On these steps in the 
/-/ plot the neuron generates n action potentials during m 
cycles of the periodic drive for a range of values of / (or in- 
trinsic frequency). Noise induces jitter in the spike time, but 
the CV on the steps is reduced compared to the CV outside 
the steps (and steps with higher m values). The entrainment 
produces a phase difference between the firing time and the 
crest of the drive. The size of the phase difference depends on 
the intrinsic frequency of the neuron. Robustness is less for 
current values close to the edge of the step (with unfavorable 
phase differences). Therefore noise will reduce the width of 
the phase-locking step (not shown). Here we found network 
parameters for which coincidence could be maintained despite 
highly variable firing rates of the neuronal populations. 

Our work, and also a recent study |]2l|], is to a large extent 
based on the recent contributions by Wang&Buzsaki [po|]. It is 
therefore important to briefly reiterate, and spell out how our 
work extends the work of Wang & Buzsaki, and how it differs 
pi]]. We included the effect of synaptic noise. We have shown 
that for the purposes of our modeling work a Gaussian white 
noise current can adequately reproduce experimental ranges 
of CV |]4^]. Biophysically realistic amounts of noise do af- 
fect the synchronization we have studied. The noise effects 
are also different from those of current heterogeneity as was 
discussed before. The issue of noise was not addressed by 
Wang & Buzsaki. Another important difference is that previ- 
ous works [ pO| , pT| | studied only strong synchronization. Here 
we have proposed that stochastic weak synchronization under- 
lies the synchronized population oscillations in the hippocam- 
pus. For this reason our robust 40Hz population rhythms were 
obtained for different values of the coupling parameters gsyn, 
Tsyn, and the driving current /, and the system size com- 
pared to previous work . In our computational work we 
actually needed a small amount of noise to obtain weak syn- 
chronization. 

The synchronization properties of large networks may be 
of some mathematical interest. Our networks are small, and 
probably the behavior can change quantitatively when increas- 
ing the network size significantly. However, in this paper we 
only addressed the question as to whether networks of physio- 
logically realistic size and connectivity can robustly synchro- 
nize. Recent experimental work suggests that interneurons 
contact on the order of 60 other interneurons Q]. For that 
reason we only vary our network size between A^ = 10 and 
1000. 

In the introduction we mentioned recent in vivo work and 
the controversies on the functional role of synchronization. 
Our work obviously does not contribute to the understanding 
of the function of synchronization. An important question is 
what kind of synchronization can be sustained in biophysi- 
cally realistic networks. Traub and coworkers 1 12, 15| , p4Hl3| ] 
looked for physiological correlates of the gamma rhythms 
using in vitro and computational experiments. Their results 
show the crucial role of inhibition, and have provided much 
of the impetus for our work. The nervous system produces, 
for some unknown reason, periodic population activity using 
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circuitry consisting of noisy and heterogeneous neurons. Our 
results establish that it is possible for inhibitory neurons to be 
the driving force for synchronization under these conditions. 
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FIG. 1. Coherence parameter k (dashed Une, left hand scale), and average network frequency / (continuous line, right hand scale) versus 
(a) noise strength D (with <7i = 0) and (b) current heterogeneity cri (with D = 0). After a transient of 1 s, (a) time-averages were computed 
over 3 s, and (b) 2 s. In (b) each point represents the results of one independent value of the current heterogeneity. In (c,d) we plot the synaptic 
drive s{t), (c) for ai = with D = 0, 0.01, 0.04, and 0.09 from top to bottom; whereas in (d) D = with ai = 0, 0.01, 0.04, and 0.09 from 
top to bottom. The curves are offset by multiples of As = 1. We used the standard set of parameters described in Methods. 
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FIG. 2. (I) Loss of synchrony due to current heterogeneities: (a) ai = 0.02, (b) 0.07, (c) 0.08, respectively. (II) Loss of synchrony due 
to synaptic noise: D — 0.01 (a), 0.04 (b), 0.09 (c), respectively. In 1 and 13 we plot the autocorrelation of Xi (spike train of neuron 1); In 

2 — 10 we plot the cross correlation of Xi with X2, ■■, Xio, respectively; In 11 we show the cross correlations between Xi and Xj, averaged 
over all pairs In 12 we plot the autocorrelation of Xi, averaged over all neurons, and in 14 the auto correlation of the total spiking rate 
X; In 15 we show the rastergrams of the network, i.e. neuron munber versus spiking time. The time-scale bar, shown in Icl2, applies to the 
curves 1 — 14, (a-c) in I and II. The y-axis is in arbitrary units, and the same scale is used for the curves in 1 — 11, and 12 — 14, except for the 
ciu^'es Ila 12,13 and lib 12,13, which are rescaled by a factor of 10. After a transient of 1 s, the time averages are computed over 2 s (I), and 

3 s (II). We used the standard set of parameters described in Methods. 
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The noise strengths are (c) D = 0.01, (d) 0.04, and (e) 0.09, respectively. In all graphs / = 1.0 except in case (4), with / = 0.61. The 
time-scale bar for the ISIH is shown in c 1, the y-scale is arbitrary, but it is the same for all curves with the same £>-value. The averaging times 
are 10 s (network), and 200 s (isolated neuron). 
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FIG. 5. (a) Coherence parameter k, (b) the coefficient of variation CV for a network, and (c) CV for an isolated neuron, versus the applied 
current / for different values of D = 0.0 (continuous line), 0.004 (dotted line), 0.04 (dashed line), and 0.2 (dot-dashed line). ISIHs for a 

network of neurons (1,3), and an isolated neuron (2,4) for different values of the current (d) / = 0.8, (e) 0.96, and (f) 1.14. We used D = 0.0 
(1,2) and 0.004 (3,4). The bar time-scale for the ISIH is shown in d2. We applied an excitatory drive with reversal potential E = mV, 
strength Qsyn = 0.05, and frequency 40 Hz. The isolated neuron has an autosynaptic feedback of Qsyn = 0.1. After a transient of 1 s (0.5 s), 
time-averages were calculated over 10 s (50 s) for the network (isolated neuron). 
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FIG. 6. (a) Network frequency, /, and (b) synchronization parameter, k, plotted versus noise strength D, for different values of the synaptic 
coupling and applied current {gsyn, /)=(0.02, 0.6955), (0.1, 1.0), (0.3, 1.625), and (0.5, 2.15), labeled by (l)-(4), respectively. The value of 
the applied current is chosen such that the neuron network will fire at / ~ 39.05 Hz at D = 0. In (c) we plot k versus Qsyn for D = 0.02 at 
the aforementioned current values. Time averages are computed over 10 s after a transient of 1 s. 
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FIG. 7. Rastergrams comparison, (a) Qayn = 0.1, / = 1.0, and (b) Qsyn = 0.5, / = 2.15, with noise strength D = 0.002. We also 
compare the ISIHs for different values of the synaptic coupling gsyn = 0.5 (1), 0.2 (2), and 0.05 (3), and for different values of the noise 
strength D equal to (c) 0.014, (d) 0.034, (e) 0.07 (e), and (f) 0.17. The values of the applied currents are / = 2.15, 1.331, and 0.8145 for (1), 
(2), and (3), respectively. Time-scale bar is shown in fl, the y-scale is arbitrary but the same for all curves in (c)-(f). For clarity the top of cl 
and dl is cut off. An initial transient of 10 s in (a) and (b) was discarded. After a transient of 1 s, time-averages were taken over 10 s (c-f). 
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FIG. 8. We plot (a) network frequency f„, (b) single neuron firing rate /, (c) cluster size A^o normalized by system size N = 100, (d) k 
, (e) coefficient of variation CV of ISI, and (f) standard deviation of the network period t„ versus synaptic strength Qsyn, for different values 
of _D = 0.0 (dot-dashed line), D — 0.008 (dotted) and D = 0.2 (solid line). gsy„ is varied from 0.05 to 2.5 in increments of 0.05. For 
D = 0.008 (D = 0.2) the first 6 (14) points have been removed according to the criterium given in the Methods section. In (b) we have added 
the firing rate of a single neuron with autosynaptic feedback (fat solid line), and the network for 13 = (open circles). Other parameters are 
CT/ = 0, / = 2.0, Tsyn = 20. A 500 ms transient is discarded, and averages are taken over 2000 ms. To smooth (d-f) we have performed 
running averages over five points, the original points are denoted by dots. 
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FIG. 9. We plot (a)-(e) as in the previous figure, (f) number Ns of active neurons for different values of -D = 0.0 (solid), D = 0.008 
(dot-dashed). Other parameters are ai = 0.02, / = 2.0, Tsy„ = 20. The smoothed (d-f) figures were obtained as in the previous figure. 
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FIG. 10. We plot (a) network frequency f„, (b) single neuron firing rate /, (c) coefficient of variation CV of ISI, (d) k , (e) cluster size 
Nc normalized by system size N, (f) number of active neurons Ns normalized by N, (g) coefficient of variation of Nc, (h) standard deviation 
of the network period r„ versus system size A'^ (A'^ = 10, 20, 50, 100, 200, 500, 1000). We show results for ten different realizations of the 
current distribution (open circles), and the their average (solid lines). Other parameters are ai = 0.10, D = 0.2, I = 5.0, gsyn = 1-0, and 
Tsyn = 20. A 500 ms transient is discarded, and averages are over at least 2000 ms. 
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FIG. 1 1. We plot in (a,c) the interspike intervals (ISI) of the network neurons sorted from the lowest to highest value, and (b,d) the ISIHs. 
In (a,c) we used different values of the noise-strength, from top to bottom D = 0.0, 0.04, 0.08, 0.36, and 0.56, with ai = 0. In (c,d) from 

top to bottom o-j = (1), 0.173 (2),0.346 (3), 0.520 (4), and 0.693 (5), with D = 0.2. After a transient of 500 ms, the time-average is 
computed over 2 s. We plot the instantaneous firing rate as a function of time for (e) D = 0, ai = 0, (f) D = 0.56, ai = 0, and (g) D = 0.2, 
ai = 0.35. The scale bars are (g) 500 imp/s, (h,i) 50 imp/s. We used / = 3.5, Tsyn = 20, N = 1000, and Qsyn = 2.0. 
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FIG. 12. We plot in (a,d) the network frequency /« (solid line, left hand scale) and single neuron firing rate / (dashed line, right hand 
scale); (b,e) coherence k (solid line) and the CV (dashed line); (c,f) cluster size Nc and number of active neurons Ng divided by JV; as function 

of (a,b,c) noise variance D, and (d,e,f) current heterogeneity ai . We used the following parameters I = 3.5, Tsyn = 20, Qsyn = 0, and 
N = 1000, for (a,b,c) with ai = and for (d,e,f) the same parameters with D = 0.2. 
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FIG. 13. We plot in (a,b) the interspike intervals (ISI) for the 100 neurons sorted from the lowest to highest values, and in (c,d) the ISlHs. 
In (a,c) we used different values for the noise-strength, from top to bottom D = 0.04, 0.24, 0.84, 4.2, and 8.2, with ai = 0. In (b,d) from top 
to bottom aj = (1), 0.04 (2),0.08 (3), 0.12 (4), and 0.16 (5), with D = 1.04. A sinusoidal current with a period of 25 ms and an amplitude 
1.2 iiA/cm^ is added to the soma. The time-scale bar is shown in dl, the y-scale is arbitrary but the same in (c), and (d). After a transient of 
2 s, the time-average is computed over 10 s. (e,f) Instantaneous spiking rate versus time for (e) crj = 0, and (f) crj = 0.10, both at D = 1.04. 
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